﻿#'
#' @TODO ## 差异分析
#' @title 差异分析
#' @description 利用R包`limma`进行差异分析。输入表达谱与分组数据，进行差异分析。.
#' @param exp 表达谱数据
#' @param group_infor 分组信息, 分组信息列名以`group`做列名,样本所在列列名为`sample`
#'  > 注： group中默认第一会作为`control`，另一分组作为`case`。
#' 例如 tumor与normal，字母n在t前面，所以第一顺位是normal，以normal 为对照，观察tumor的上调与下调
#' @param output_dir 输出文件目录
#' @param var_name 用于文件命名
#' @param levels_manual 指定差异分析中的levels，如果是NULL，则默认group_infor中第一顺位的分组为对照组，第二个为实验组。
#' 否则指定一个顺序，是group_infor中的组别名称。eg：levels_manual = c('Normal','Tumor')
#'
#' @return  a *list*，第一个元素全部degs结果，第二个为表达谱，是按分组文件筛选后的表达谱
#' @export
#'
#' @author *WYK*
#'
DEGs <- function(exp = NULL, group_infor = NULL, output_dir = "./", var_name = NULL, levels_manual = NULL) {
  if (is.null(levels_manual)) {
    levels_manual <- c(sort(unique(group_infor$group))[1], sort(unique(group_infor$group))[2])
  }

  if (length(var_name) == 0) {
    var_name <- paste0(sample(c(letters, LETTERS), 5), collapse = "")
  }

  suppressMessages(library(magrittr))

  library(tidyverse)
  group_infor <- na.omit(group_infor)
  rownames(group_infor) <- NULL

  if(!'group' %>% is_in(colnames(group_infor))){
    warning("缺少分组信息列，或者 分组信息所在列列名不是group，程序跳出。")
    return()
  }

  common_sample <- intersect(group_infor$sample, colnames(exp))
  group_infor <- group_infor %>%
    column_to_rownames("sample") %>%
    .[common_sample, , drop = F] %>%
    rownames_to_column("sample")
  exp_1 <- exp[, common_sample]

  colnames(exp_1) == group_infor$sample

  library(limma)
  eset <- as.matrix(exp_1)

  group <- factor(group_infor$group, levels = levels_manual, ordered = F)
  design <- model.matrix(~group)
  colnames(design) <- levels(group)

  fit <- lmFit(eset, design)
  fit2 <- eBayes(fit)
  dif <- topTable(fit2, coef = 2, n = nrow(fit2))
  # head(dif)
  degs <- dif %>%
    as.data.frame() %>%
    rownames_to_column(var = "gene") %>%
    dplyr::rename(log2FC = logFC, pvalue = P.Value, adjusted_pvalue = adj.P.Val) %>%
    dplyr::select(gene, log2FC, pvalue, adjusted_pvalue)

  if (!dir.exists(sprintf("%s/output/", output_dir))) {
    dir.create(sprintf("%s/output/", output_dir), recursive = T)
  }

  degs %>% write_tsv(
    x = ., file = sprintf("%s/output/deg_res_%s.txt", output_dir, var_name),
    quote = "none",
  )

  return(degs)
}


#' @title  ## 按阈值筛选差异分析结果
#' @TODO  按阈值筛选差异分析结果
#' @description 按设定阈值对差异分析结果进行分析。差异分析全部结果与各项阈值.
#' @param degs gene、log2FC等数据
#' @param log2FC_value log2FC绝对值，默认1，即两倍差异
#' @param DEG_adjusted_pvalue 矫正后p值筛选结果，默认0.05
#' @param DEG_pvalue 如果`DEG_adjusted_pvalue`不支持后续分析，可以使用p值筛选结果。
#'   > 注：如果`DEG_pvalue != NULL`，在程序运行中，其优先级高于`DEG_adjusted_pvalue`，
#'   >     即`DEG_pvalue`与`DEG_adjusted_pvalue`都存在的情况下，会优先使用`DEG_pvalue`
#' @param intrest_gene 兴趣基因集，字符串向量,如果为NULL，则不对差异结果进行筛选
#' @param exp_1 上一步DEGs，筛选后的表达谱数据
#'
#' @export
#' @author *WYK*
#'
###
model_data <- function(degs = NULL, log2FC_value = 1, DEG_adjusted_pvalue = 0.05,
                       DEG_pvalue = NULL, intrest_gene = NULL) {
  if (is.null(DEG_pvalue)) {
    deg1 <- degs %>%
      filter(adjusted_pvalue < DEG_adjusted_pvalue) %>%
      filter(log2FC > log2FC_value | log2FC < (-log2FC_value))
  } else {
    deg1 <- degs %>%
      filter(pvalue < DEG_pvalue) %>%
      filter(log2FC > log2FC_value | log2FC < (-log2FC_value))
  }

  up <- deg1 %>%
    filter(log2FC > log2FC_value) %>%
    nrow(.)
  down <- deg1 %>%
    filter(log2FC < -log2FC_value) %>%
    nrow(.)


  if (is.null(DEG_pvalue)) {
    msg_tmp <- sprintf(
      "There're %s  DEgenes with log2FC = %s, adj.p = %s. Up:%s , Down:%s .", nrow(deg1),
      log2FC_value, DEG_adjusted_pvalue, up, down
    )
    print(msg_tmp)
  } else {
    msg_tmp <- sprintf(
      "There're %s  DEgenes with log2FC = %s, p = %s. Up:%s , Down:%s .", nrow(deg1),
      log2FC_value, DEG_pvalue, up, down
    )
    print(msg_tmp)
  }


  if (length(intrest_gene != 0)) {
    deg2 <- deg1[which(deg1$gene %in% intrest_gene), ]
    if (nrow(deg2) == 0) {
      sprintf(
        "In \"%s\"'s setting, 0 gene mapped in gene list.",
        msg_tmp
      )
    } else {
      msg_tmp_1 <- sprintf("In \"%s\"'s setting, %s gene mapped in gene list.", msg_tmp, nrow(deg2))
      message(msg_tmp_1)
    }
  } else {
    deg2 <- deg1
  }

  # tmp <- list(deg2 = deg2, exp_2 = exp_2)

  return(deg2)
}


#' @TODO  Univariate_Cox
#' @title   ## Univariate_Cox
#'
#' @description  使用二分类变量，计算HR，与其他所需p值以及95CI.
#' @param by_group 是否使用二分类变量进行单因素cox回归分析，默认F
#' @param exp 格式要求：表达谱数据，基因在行样本在列的data.frame
#' @param clinical 临床信息数据 必须`sample`做样本列名，`time`做生存时间列名，`status`做生存状态列名
#' @param pval_univ_cox  单因素结果，p值阈值，默认0.05。如果是1，则保留全部结果
#' @param best 布尔值，是否使用最优分组进行分析
#' @return  a *tibble*
#' @export
#' @usage
#'    HR_data <- univariate_Cox(exp = ssgsea_score, clinical = sur_data, pval_univ_cox = 0.05，by_group = F)
#'    #单因素cox回归计算HR
#' @author *WYK*
#'
univariate_Cox <- function(exp = NULL, clinical = NULL, pval_univ_cox = 1, by_group = F, best = FALSE) {
  library(tidyverse)
  library(survival)
  library(survminer)
  inner_join_q <- purrr::quietly(inner_join)

  clinical <- clinical %>% filter(time > 0) %>% distinct()

  clinical <- clinical[na.omit(match(colnames(exp), clinical$sample)), ]
  exp <- exp[, na.omit(match(clinical$sample, colnames(exp)))]

  i <- which(rowSums(exp) == 0)
  if(length(i)> 0){
    cli::cli_alert_warning("有基因表达全为0，自动删除, {rownames(exp)[i]}")
    exp <- exp[-i,]
  }
  
  exp <- as.data.frame(t(exp))
  gene_name <- colnames(exp)

  if (isFALSE(best) && by_group) {
    exp %<>% mutate(across(everything(), ~ ifelse(.x > median(.x), 1, 0)))
  }

  i <- 1
  df1 <- tibble()
  df2 <- tibble()
  result <- data.frame()

  if (isTRUE(best)) {
    for (i in gene_name) {
      # cat(i)
      if (length(which(exp[, i] != 0)) < as.integer(length(exp[, i]) / 6)) {
        exp <- exp[, -which(colnames(exp) == i)]
        gene_name <- gene_name[-which(gene_name == i)]
        cli::cli_alert_warning("remove gene {i}, because if expression in too small samples.")
      }
    }
  }

  result2 <- sapply(gene_name, function(gene_name) {
    exp %>%
      rownames_to_column("sample") %>%
      inner_join_q(clinical) %>% .$result -> d

    if (isTRUE(best)) {
      sur.cut <- surv_cutpoint(
        data = d, time = "time", event = "status",
        variables = gene_name
      )
      cut <- summary(sur.cut)$cutpoint

      d[[gene_name]] <- ifelse(d[[gene_name]] > cut, 1, 0)
  
    }

    result <- coxph(
      formula = Surv(time, status) ~ get0(gene_name), data = d
    )

    result <- summary(result)
    coef <- result$coef[1]
    Hazard_Ratio <- result$coef[2]
    p.value <- result$wald["pvalue"]
    lower_.95 <- result$conf.int[, "lower .95"]
    upper_.95 <- result$conf.int[, "upper .95"]
    logrank_pvalue <- result$sctest["pvalue"]
    wald_pvalue <- result$waldtest["pvalue"]
    Likelihood_pvalue <- result$logtest["pvalue"]

    df1 <- as_tibble(
      cbind(gene_name, coef, p.value, Hazard_Ratio, lower_.95, upper_.95, logrank_pvalue, wald_pvalue, Likelihood_pvalue)
    )
    df2 <- bind_rows(df2, df1)
    return(df2)
  })

  result2 <- as.data.frame(t(result2))
  result2[, 2:ncol(result2)] <- apply(result2 %>% dplyr::select(2:ncol(result2)), 2, as.numeric)

  colnames(result2)[1] <- "gene"
  result2 <- result2 %>%
    filter(p.value < pval_univ_cox)
  result2$gene <- unlist(result2$gene)

  result2 <- result2 %>%
    mutate(HR = paste0(
      sprintf("%.2f", Hazard_Ratio),
      "(",
      sprintf("%.2f", lower_.95),
      "-",
      sprintf("%.2f", upper_.95),
      ")"
    )) %>%
    as_tibble()

  return(result2)
}


#' @title ## 计算lasso结果
#' @TODO: 计算lasso结果
#' @description   调用`compute_lasso_coef`子函数，计算score，绘制lasso回归系数图与coef条形图.
#' 
#' lasso回归的结果解读：LASSO 回归的特点是在拟合广义线性模型的同时进行变量筛选和复杂度调整。
#' 其中，复杂度调整的程度由参数 lambda来控制， lambda越大，对变量较多的线性模型的惩罚力度就越大，从而最终获得一个变量较少的模型。
#' 首先我们看第一个结果图，即Lasso筛选变量过程图，自变量运动轨迹图
#'
#' - 图中的每一条曲线代表了每一个自变量系数的变化轨迹。
#'   纵坐标是系数的值，下横坐标是 log(λ) ，上横坐标是此时模型中非零系数的个数。
#'   黑线代表的自变量 1 在 λ值很大时就有非零的系数，然后随着 λ 值变小不断变大。
#' 
#' - 交叉验证图,
#'   交叉验证，对于每一个 λ 值，在红点所示目标参量的均值左右，我们可以得到一个目标参量的置信区间。
#'   两条虚线分别指示了两个特殊的 λ值
#'   lambda.min 是指在所有的 λ 值中，得到最小目标参量均值的那一个。 
#'   而 lambda.1se 是指在 lambda.min 一个方差范围内得到最简单模型的那一个 λ 值。
#'   因为 λ 值到达一定大小之后，继续增加模型自变量个数即缩小 λ 值，并不能很显著的提高模型性能， lambda.1se 给出的就是一个具备优良性能但是自变量个数最少的模型。
#' @param clinical  clinical_data生存信息，表头**必须**包含time与status，
#' @param exp  基因在行，样本在列的表达谱，
#' @param unicox_res_gene  单因素回归后，筛选得到的基因
#' @param seed  默认种子为1110
#' @param saveplot saveplot是否在当前文件创建lasso文件夹保存图片
#' @param output_dir  lasso结果图片输出目录
#' @export
#' @return `list`
#' @return
#'  - 返回一个list，第一个元素为geneSymbol与coef组成的data.frame，第二个元素为fit，可用于plot(),第三个元素为cvfit，可用于plot()。
#'  - 图片描述：A: 每个自变量的变化轨迹，横轴代表自变量 lambda 的log 值，纵轴表示自变量的系数；B: 每个 lambda下的置信区间；C: 特征基因的回归系数 From llk
#' @author *WYK*
#'
#'
lasso <- function(exp = NULL, clinical = NULL, unicox_res_gene = NULL, seed = 1110, saveplot = F, output_dir = "./",
                  var_name = NULL, height_used = 4.75) {
  if (length(var_name) == 0) {
    var_name <- deparse(substitute(exp))
  }

  lasso_res <- compute_lasso_coef(
    clinical = clinical, exp = exp,
    unicox_res_gene = unicox_res_gene, seed = seed
  )

  temp_coef <- lasso_res[[1]] %>%
    mutate(
      group = ifelse(.$coef > 0, "1", "2"),
      min_coef = min(coef), max_coef = max(coef)
    )

  fit <- lasso_res[[2]]
  cvfit <- lasso_res[[3]]

  # x <- coef(fit)
  # tmp <- as.data.frame(as.matrix(x))
  # tmp$coef <- row.names(tmp)
  # tmp <- reshape2::melt(tmp, id = "coef")
  # tmp$variable <- as.numeric(gsub("s", "", tmp$variable))
  # tmp$coef <- gsub('_','-',tmp$coef)
  # tmp$lambda <- fit$lambda[tmp$variable+1] # extract the lambda values
  # tmp$norm <- apply(abs(x[-1,]), 2, sum)[tmp$variable+1] # compute L1 norm

  # plot(fit,xvar = 'lambda',label = F)
  # plot(cvfit)

  data_for_fit <- fit %>%
    coef() %>%
    as.matrix() %>%
    as.data.frame() %>%
    rownames_to_column(var = "gene") %>%
    pivot_longer(cols = -gene, names_to = "variable", values_to = "value") %>%
    mutate(variable = as.numeric(gsub("s", "", variable))) %>%
    mutate(gene = str_replace_all(gene, "-", "_")) %>%
    mutate(lambda = fit$lambda[.$variable + 1])

  library(ggsci)

  fit_plot <- ggplot(data_for_fit, aes(log(lambda), value, color = gene)) +
    geom_vline(xintercept = log(cvfit$lambda.min), size = 0.8, color = "grey60", alpha = 0.8, linetype = 2) +
    geom_line(size = 1) +
    xlab("log(Lambda)") +
    ylab("Coefficients") +
    egg::theme_article(base_size = 12) +
    theme(plot.background = element_rect(fill = "white")) +
    scale_color_manual(values = c(
      pal_npg()(10),
      pal_d3()(10),
      pal_jco()(7),
      pal_lancet()(9),
      pal_nejm()(8),
      unique(unlist(lapply(
        as.list(1:2600),
        function(i) paste0(sample(0:255, size = 3, replace = TRUE), collapse = " ")
      )))
    )) +
    scale_x_continuous(expand = c(0.01, 0.01)) +
    scale_y_continuous(expand = c(0.01, 0.01)) +
    theme(
      panel.grid = element_blank(),
      axis.title = element_text(size = 15, color = "black"),
      axis.text = element_text(size = 12, color = "black"),
      legend.title = element_blank(),
      legend.text = element_text(size = 12, color = "black"),
      legend.position = "none"
    ) +
    annotate("text",
      x = max(log(data_for_fit$lambda)) * 1,
      y = min(data_for_fit$value) * .85,
      label = sprintf("Optimal Lambda = %.4f", cvfit$lambda.min),
      color = "black", hjust = "right", size = 4.35
    )

  data_for_cvfit <- tibble(
    lambda = cvfit[["lambda"]], cvm = cvfit[["cvm"]], cvsd = cvfit[["cvsd"]],
    cvup = cvfit[["cvup"]], cvlo = cvfit[["cvlo"]], nozezo = cvfit[["nzero"]]
  ) %>%
    mutate(log_lambda = log(lambda))

# sec_xlabel <- approx(x=index,y=df,xout=xbreaks,rule=2,method="constant",f=1)$y

  cvfit_plot <- ggplot(data_for_cvfit, aes(log_lambda, cvm)) +
    geom_errorbar(aes(x = log_lambda, ymin = cvlo, ymax = cvup), width = 0.05, size = 0.6, alpha = 0.5) +
    geom_vline(
      xintercept = data_for_cvfit$log_lambda[which.min(data_for_cvfit$cvm)], size = 0.8,
      color = "grey60", alpha = 0.8, linetype = 2
    ) +
    geom_vline(
      xintercept = log(cvfit$lambda.1se), size = 0.8,
      color = "grey60", alpha = 0.8, linetype = 2
    ) +
    geom_point(size = 2, color = c("#c41e1e")) +
    xlab("log(Lambda)") +
    ylab("Partial Likelihood Deviance") +
    egg::theme_article(base_size = 12) +
    theme(plot.background = element_rect(fill = "white")) +
    scale_color_manual(values = c(
      pal_npg()(10),
      pal_d3()(10),
      pal_jco()(7),
      pal_lancet()(9),
      pal_nejm()(8)
    )) +
    scale_x_continuous(
      expand = c(0.02, 0.02)
    ) +
    scale_y_continuous(expand = c(0.02, 0.02)) +
    theme(
      panel.grid = element_blank(),
      axis.title = element_text(size = 15, color = "black"),
      axis.text = element_text(size = 12, color = "black"),
      legend.title = element_blank(),
      legend.text = element_text(size = 12, color = "black"),
      legend.position = "none"
    ) +
    annotate("text",
      x = log(cvfit$lambda.min),
      y = max(data_for_cvfit$cvup) * 1.02,
      label = sprintf("var = %s", nrow(lasso_res[[1]])),
      color = "black", hjust = "center", size = 4.35
    )

  # annotate('text',x = data_for_cvfit$log_lambda[which.min(data_for_cvfit$cvm)],
  #                     y = max(data_for_cvfit$cvup)*1.004,
  #          label = data_for_cvfit$nozezo[which.min(data_for_cvfit$cvm)],size = 5)

  # annotate('text',x = max(log(data_for_fit$lambda))*1,
  #          y=max(data_for_cvfit$cvup)*1.004,
  #          label=sprintf('Optimal Lambda = %.4f',cvfit$lambda.min),
  #          color='black',hjust = 'right',size = 4.35)

  coef_bar_p <- ggplot(
    data = temp_coef,
    mapping = aes(x = coef, y = factor(reorder(symbol, coef)), fill = factor(group))
  ) +
    geom_bar(stat = "identity", color = "black", width = .8, size = .6) +
    theme_pubr() +
    guides(fill = "none") +
    ylab("Gene") +
    xlab("Coefficients") +
    # theme(
    #   axis.ticks.y = element_blank(),
    #   axis.line.y = element_blank()
    # ) +
    # scale_x_continuous(
    #   expand = c(0, 0),
    #   breaks = seq(round(temp_coef$min_coef[1] - 0.1, 1), round(temp_coef$max_coef[1] + 0.1, 1), 0.4)
    # ) +   
    scale_fill_manual(values = c("1" = "#3171a5", "2" = "#FF7F00"))
  # ggfittext::geom_bar_text(data = temp_coef %>% filter(group == '1'),
  #                          aes(x = 0,label = sprintf('%.4f',coef)),outside = T,place = 'left',reflow = T)+
  # ggfittext::geom_bar_text(data = temp_coef %>% filter(group == '2'),
  #                          aes(x = 0,label = sprintf('%.4f',coef)),outside = T,place = 'right')


  lasso_res[[4]] <- coef_bar_p
  p <- cowplot::plot_grid(
    fit_plot + theme(plot.margin = unit(c(1.25, 0.1, 1.25, 0.2), "cm")),
    cvfit_plot + theme(plot.margin = unit(c(1.25, 0.1, 1.25, 0.2), "cm")),
    coef_bar_p + theme(plot.margin = unit(c(.85, 0.35, 0.25, 0.1), "cm")), # 调整画布大小,
    labels = "AUTO",
    label_size = 20, nrow = 1, rel_widths = c(.97, .97, 1)
  )

  lasso_res[[5]] <- p

  # print(p)


  if (saveplot) {
    if (!dir.exists(sprintf("%s/output/lasso_res/", output_dir))) {
      dir.create(sprintf("%s/output/lasso_res/", output_dir), recursive = T)
    } else {
      (
        print("Dir is ready.")
      )
    }

    dir_now <- paste0(output_dir, "/output/lasso_res/")

    lasso_res[[1]] %>% write_delim(x = ., file = sprintf("%slasso_coef_%s.txt", dir_now, var_name), delim = "\t", quote = "none")

    cowplot::ggsave2(
      filename = sprintf("%slasso_%s.pdf", dir_now, var_name),
      plot = p, width = 12, height = height_used
    )
    cowplot::ggsave2(
      filename = sprintf("%slasso_%s.tiff", dir_now, var_name),
      plot = p, width = 12, height = height_used, dpi = 300
    )
    cowplot::ggsave2(
      filename = sprintf("%slasso_%s_72dpi.tiff", dir_now, var_name),
      plot = p, width = 12, height = height_used, dpi = 72
    )
  }

  names(lasso_res) <- c("coef", "fit", "cv.fit", "Fig_coef_bar", "Fig_all")

  return(lasso_res)
}
# lasso子函数，主要用于计算coef结果
#' @export
compute_lasso_coef <- function(exp = NULL, clinical = NULL, unicox_res_gene = NULL, seed = 1110) {
  library(cowplot)
  library(tidyverse)
  library(survival)
  library(ggpubr)
  var_name <- deparse(substitute(exp))

  unicox_exp <- exp[unicox_res_gene, ]
  clinical <- clinical[na.omit(match(colnames(unicox_exp), clinical$sample)), ]
  unicox_exp <- unicox_exp[, na.omit(match(clinical$sample, colnames(unicox_exp)))]
  # dim(unicox_exp)
  # dim(clinical)

  lasso_res <- vector('list',3L)
  library(glmnet)
  y <- data.matrix(Surv(clinical$time, clinical$status))
  x <- unicox_exp %>%
    t() %>%
    as.data.frame() %>%
    as.matrix()

  set.seed(seed)
  fit <- glmnet(x, y, family = "cox", alpha = 1)
  lasso_res[[2]] <- fit
  # plot(fit)

  set.seed(seed)
  cvfit <- cv.glmnet(x, y, family = "cox")
  lasso_res[[3]] <- cvfit
  # plot(cvfit)

  tmp <- coef(object = cvfit, s = "lambda.min")

  Signature_Coef_min_os <- tmp %>%
    as.matrix() %>%
    as.data.frame() %>%
    rownames_to_column("gene") %>%
    dplyr::rename(coef = 2) %>%
    filter(coef != 0)

  colnames(Signature_Coef_min_os) <- c("symbol", "coef")

  message(paste0("Univ_cox has ", nrow(unicox_exp), " genes"),'\n')
  message(paste0("Lasso-cox gets ", dim(Signature_Coef_min_os)[1], " genes"))

  lasso_res[[1]] <- Signature_Coef_min_os

  return(lasso_res)
}


#' @TODO 得到KM曲线与ROC曲线
#' @title ## 数据集KM曲线与ROC曲线
#' @description 依据`lasso`模型计算riskscore、并高低分组集成数据框,绘制KM、ROC曲线
#' 主函数：KM_ROC_curve，包含三个子函数:
#' - `compute_score`用于计算riskscore,
#' - `get_group`用于计算得到含有score以及高低分组信息的临床信息表
#' - `KM_ROC`利用上一个子函数得到的结果包含高低风险分组信息的矩阵绘制KM曲线与ROC曲线。
#'
#' 会在output_dir路径下生成一个output文件夹，output文件夹内，会以exp变量名命名KM-ROC曲线文件夹名字.
#'
#' @param model_coef  包含系数的文件，第一列为基因名，第二列为具体coef
#' @param exp 整理后的表达数据
#' @param clinical 整理后的临床信息，需要以`sample`、`time`、`status`为列名，对应样本名、生存时间、生存状态
#' @param surtime_unit 生存信息中的时间单位，依据临床信息中的time单位，可选1、12或者365，对应时间单位是年、月、日
#' @param saveplot 是否保存图片
#' @param ROC_time_break ROC曲线时间断点，默认1,3,5年
#' @param best_cut 逻辑值，是否是否使用最优cutoff来绘制KM曲线
#' @param do_ROC_CI 逻辑值，是否是在图中给出ROC曲线的的95%CI
#' @param manual_cutoff 数值，分组阈值
#' @param savetiff 逻辑值，在`saveplot = T`的条件下，是否生成与保存tiff文件
#' @export
#' @return *list*
#'  - 元素1：ggsurvplot对象KM曲线
#'  - 元素2：gg对象ROC曲线
#'  - 元素3：结果data.frame
#'  - 元素4：p_auc_res p值，AUC等信息
#' @return reslut <- KM_ROC_curve(……)，自定义名字，可选`saveplot = T`，保存KM与roc曲线
#' @author *WYK*
#'
KM_ROC_curve <- function(model_coef = NULL, exp = NULL, clinical = NULL,
                         surtime_unit = c(12, 365), saveplot = T, output_dir = "./", var_name = NULL,
                         ROC_time_break = c(1, 3, 5), best_cut = F, tri = F,
                         manual_cutoff = NULL, do_ROC_CI = F, savetiff = F) {
  exp <- as.data.frame(exp)
  clinical <- as.data.frame(clinical)

  riskscore_res <- compute_score(model_coef = model_coef, exp = exp, clinical = clinical)
  score_res_data <- get_group(riskscore_res = riskscore_res, clinical = clinical, exp = exp)
  KM_ROC_curve <- KM_ROC(
    score_res = score_res_data, surtime_unit = surtime_unit,
    ROC_time_break = ROC_time_break, best_cut = best_cut, manual_cutoff = manual_cutoff,tri = tri,
    do_ROC_C = do_ROC_CI
  )
  KM_ROC_curve[[3]] <- KM_ROC_curve[[3]] %>% as.data.frame()

  if (length(var_name) == 0) {
    var_name <- paste0(sample(c(letters, LETTERS), 5), collapse = "")
  }

  if (saveplot) { # & !is.na(KM_ROC_curve[[1]])

    if (!dir.exists(sprintf("%s/output/KM_ROC_curve_%s", output_dir, var_name))) {
      dir.create(sprintf("%s/output/KM_ROC_curve_%s", output_dir, var_name), recursive = T)
    } else {
      print(sprintf("Dir '%s/output/KM_ROC_curve_%s' is existed.", output_dir, var_name))
    }

    plot_dir <- sprintf("%s/output/KM_ROC_curve_%s", output_dir, var_name)

    p <- ggarrange(KM_ROC_curve[[1]]$plot + labs(x = "") + theme(plot.margin = unit(c(.2, .2, -.165, .2), "cm")),
      KM_ROC_curve[[1]]$table + theme(plot.margin = unit(c(-.15, .2, .2, .2), "cm")),
      ncol = 1, align = "v", heights = c(.73, .27)
    )

    write_tsv(x = KM_ROC_curve[[3]], file = sprintf("%s/table_%s_score_res.txt", plot_dir, var_name), quote = "none")

    cowplot::ggsave2(
      filename = sprintf("%s/figure_%s_KM.pdf", plot_dir, var_name),
      plot = p, width = 5, height = 5.5
    )

    cowplot::ggsave2(
      filename = sprintf("%s/figure_%s_ROC.pdf", plot_dir, var_name),
      plot = KM_ROC_curve[[2]], width = 5.5, height = 5.2
    )


    p2 <- cowplot::plot_grid(p, KM_ROC_curve[[2]],
      nrow = 1, labels = "AUTO", rel_widths = c(0.95, 1.05), label_size = 20, scale = c(1, .97)
    )

    cowplot::ggsave2(
      filename = sprintf("%s/figure_%s_KM_ROC.pdf", plot_dir, var_name),
      plot = p2, width = 10, height = 5.2
    )

    if (savetiff) {
      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_KM.tiff", plot_dir, var_name),
        plot = p, width = 5.2, height = 5.5, dpi = 300
      )
      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_KM_dpi72.tiff", plot_dir, var_name),
        plot = p, width = 5.2, height = 5.5, dpi = 72
      )

      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_ROC.tiff", plot_dir, var_name),
        plot = KM_ROC_curve[[2]], width = 5.5, height = 5.2, dpi = 300
      )
      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_ROC_dpi72.tiff", plot_dir, var_name),
        plot = KM_ROC_curve[[2]], width = 5.5, height = 5.2, dpi = 72
      )

      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_KM_ROC.tiff", plot_dir, var_name),
        plot = p2, width = 10, height = 5.2, dpi = 300
      )
      cowplot::ggsave2(
        filename = sprintf("%s/figure_%s_KM_ROC_dpi72.tiff", plot_dir, var_name),
        plot = p2, width = 10, height = 5.2, dpi = 72
      )
    }
  }

  names(KM_ROC_curve) <- c("KM_p", "ROC_p", "ScoreInfor", "p_auc")

  return(KM_ROC_curve)
}
# 函数KM_ROC_curve的子函数，封装score计算函数
#' @export
compute_score <- function(exp = NULL, model_coef = NULL, clinical = NULL) {
  library(tidyverse)

  if (!exists("var_name")) {
    var_name <- deparse(substitute(exp))
  }

  if (length(clinical) != 0) {
    clinical <- clinical[na.omit(match(colnames(exp), clinical[, "sample"])), ] %>%
      dplyr::filter(time > 0 & time != "NA")
    exp <- exp[, na.omit(match(clinical[, "sample"], colnames(exp)))]
    clinical <- clinical[na.omit(match(colnames(exp), clinical[, "sample"])), ]
  }

  exp <- exp[model_coef[, 1], ] %>%
    t() %>%
    as.data.frame()

  miss_gene <- model_coef[, 1][which(is.na(exp[1, ]))]
  miss_gene_symbol <- paste0(miss_gene, collapse = " ")

  exp[is.na(exp)] <- 0

  # print(dim(exp))
  # print(dim(model_coef[, 2, drop = F]))

  riskscore_res <- as.matrix(exp) %*% as.matrix(model_coef[, 2, drop = F]) %>%
    as.data.frame() %>%
    dplyr::rename(riskscore = 1)

  if (length(miss_gene) == 1) {
    print(sprintf("%s is missed in %s", miss_gene_symbol, "expression"))
  } else if (length(miss_gene) > 1) {
    print(sprintf("%s are missed in %s", miss_gene_symbol, "expression"))
  } else {
    print(sprintf("All ModelGenes are mapped in %s.", "expression"))
  }

  # print(head(riskscore_res, 5))
  return(riskscore_res)
}
#' @export
# 函数KM_ROC_curve的子函数
get_group <- function(riskscore_res = NULL, exp = NULL, clinical = NULL) {
  library(tidyverse)

  if (length(exp) != 0 && length(clinical) != 0) {
    clinical <- clinical[na.omit(match(colnames(exp), clinical$sample)), ] %>%
      filter(time > 0 && time != "NA")
    exp <- exp[, na.omit(match(clinical$sample, colnames(exp)))]
    clinical <- clinical[na.omit(match(colnames(exp), clinical$sample)), ]
  } # 仅用于确认，在临床信息表中的样本，是同时具有临床信息和表达数据的样本

  score_res <- riskscore_res %>%
    rownames_to_column(var = "sample") %>%
    mutate(riskgroup = ifelse(riskscore > median(riskscore), "High", "Low")) %>%
    inner_join(clinical, .) %>%
    mutate(time = as.numeric(.$time)) %>%
    filter(time > 0) %>%
    dplyr::select(sample, status, time, riskgroup, riskscore)

  # print(head(score_res, 5))
  return(score_res)
}
#' @TODO   KM ROC曲线
#' @title  函数KM_ROC_curve的子函数
#' @param surtime_unit 时间单位，月份为12，天为365 ,年为1
#' @param ROC_time_break ROC时间断点，默认1,3,5年
#' @param score_res 包含sample，time,status，riskscore等列的得分计算结果
#' @return *list*
#'  - 元素1：KM曲线ggplot2对象
#'  - 元素2：ROC曲线ggplot2对象
#'  - 元素4：p_auc_res p值，AUC等信息
#' @export
#' @author *WYK*
#'
KM_ROC <- function(score_res = NULL, surtime_unit = c(1, 12, 365), ROC_time_break = c(1, 3, 5),
                   best_cut = F, tri = F, manual_cutoff = NULL,
                   do_ROC_CI = F) {
  library(tidyverse)
  library(survival)
  library(survminer)

  if (best_cut) {
    sur.cut <- surv_cutpoint(data = score_res, time = "time", event = "status", variables = "riskscore", minprop = 0.2)
    cut <- summary(sur.cut)$cutpoint

    score_res <- score_res %>%
      mutate(HR_group = ifelse(.$riskscore > cut, 1, 0))
  } else {
    score_res <- score_res %>%
      mutate(HR_group = ifelse(.$riskscore > median(.$riskscore), 1, 0))
  }

  if (!is.null(manual_cutoff)) {
    score_res <- score_res %>%
      mutate(HR_group = ifelse(.$riskscore > manual_cutoff, 1, 0))
  }

  if (isTRUE(tri)) {
    # score_res <- score_res %>%
    #   mutate(HR_group = ifelse(.$riskscore > manual_cutoff, 1, 0))

    cut_v <- quantile(score_res$riskscore, probs = seq(0, 1, 0.25)) %>%
      as.numeric() %>% sort() %>% 
      .[c(2, 4)]

    i1 <- which(score_res$riskscore < cut_v[1])
    i2 <- which(score_res$riskscore > cut_v[2])

    Group <- rep("Middle", time = length(score_res$riskscore))
    Group[i1] <- "Low"
    Group[i2] <- "High"

    Group <- factor(Group, levels = c("Low", "Middle", "High"))

    score_res %<>% mutate(HR_group = Group)
  }

  if (length(unique(score_res$HR_group)) != 1) {
    coxtmp <- summary(coxph(Surv(time, status) ~ HR_group, data = score_res))

    HR <- coxtmp$coefficients[2]
    logrank_pvalue <- coxtmp$sctest["pvalue"]
    lower_.95 <- coxtmp$conf.int[, "lower .95"]
    upper_.95 <- coxtmp$conf.int[, "upper .95"]
    C <- coxtmp$concordance[1]

    p_chara <- paste0(
      ifelse(logrank_pvalue < 0.001, "P < 0.001", paste0("P = ", round(logrank_pvalue, 3))),
      "\n",
      "HR = ", round(HR, 2),
      "\n95% CI = ", round(lower_.95, 2), " - ", round(upper_.95, 2)
      # "\nC-index = ", round(C, 2)
    )

    if (surtime_unit == 12) {
      chara_xlab <- "Time (Months)"
    } else if (surtime_unit == 365) {
      chara_xlab <- "Time (Days)"
    } else if (surtime_unit == 1) {
      chara_xlab <- "Time (Years)"
    } else {
      chara_xlab <- "Time"
    }


    score_res <- score_res %>%
        mutate(riskgroup = ifelse(HR_group == 1, "High", "Low"))

    if (isTRUE(tri)) {
      p_chara <- TRUE

      score_res <- score_res %>%
        mutate(riskgroup = HR_group)

      HR <- 1.1
    }

    fit <- surv_fit(Surv(time, status) ~ riskgroup, data = score_res)
    KM <- ggsurvplot(fit,
      data = score_res,
      surv.median.line = "hv",
      legend.title = "Score",
      legend.labs = c("High", "Low"),
      palette = "Set1",
      ggtheme = egg::theme_article(14), #
      pval = p_chara,
      pval.size = 4.5,
      xlab = chara_xlab,
      tables.height = 0.28,
      risk.table = T,
    )

    if (HR < 1) {
      score_res_roc <- score_res %>% mutate(status = ifelse(score_res$status == 0, 1, 0))
    } else {
      score_res_roc <- score_res
    }

    library(timeROC)
    set.seed(1110)
    ROC_data <- timeROC(
      T = score_res_roc$time,
      delta = score_res_roc$status,
      marker = score_res_roc$riskscore,
      cause = 1,
      weighting = "marginal",
      times = c(
        ROC_time_break[1] * surtime_unit,
        ROC_time_break[2] * surtime_unit,
        ROC_time_break[3] * surtime_unit
      ),
      ROC = TRUE,
      iid = F
    )

    rocFORplot <- tibble(
      year1x <- ROC_data$FP[, 1],
      year1y <- ROC_data$TP[, 1],
      year2x <- ROC_data$FP[, 2],
      year2y <- ROC_data$TP[, 2],
      year3x <- ROC_data$FP[, 3],
      year3y <- ROC_data$TP[, 3],
    )
    colnames(rocFORplot) <- c("year1x", "year1y", "year2x", "year2y", "year3x", "year3y")

    AUC_anno_1 <- sprintf("AUC at %s year = %s", ROC_time_break[1], sprintf("%.3f", ROC_data$AUC[[1]]))
    AUC_anno_2 <- sprintf("AUC at %s year = %s", ROC_time_break[2], sprintf("%.3f", ROC_data$AUC[[2]]))
    AUC_anno_3 <- sprintf("AUC at %s year = %s", ROC_time_break[3], sprintf("%.3f", ROC_data$AUC[[3]]))

    ROC_1 <- ggplot(data = rocFORplot) +
      # geom_smooth(aes(x = year1x, y = year1y), size = 1.2, color = "#E83D3F",se = F) +
      # geom_smooth(aes(x = year2x, y = year2y), size = 1.2, color = "#377EB8",se = F) +
      # geom_smooth(aes(x = year3x, y = year3y), size = 1.2, color = "#4DAF4A",se = F) +
      geom_line(aes(x = year1x, y = year1y), size = 1.2, color = "#E83D3F") +
      geom_line(aes(x = year2x, y = year2y), size = 1.2, color = "#377EB8") +
      geom_line(aes(x = year3x, y = year3y), size = 1.2, color = "#4DAF4A") +
      geom_abline(slope = 1, intercept = 0, color = "grey", size = 1, linetype = 2) +
      # plotROC::style_roc(
      #   xlab = "1-Specificity",
      #   ylab = "Sensitivity",
      #   theme = theme_test
      # ) +
      theme_test() +
      # theme(plot.background = element_rect(fill = "white")) +
      annotate(geom = "line", x = c(.43, 0.49), y = .17, colour = "#E83D3F", size = 1.2) +
      annotate("text", x = 0.5, y = .17, size = 4.35, label = AUC_anno_1, color = "black", hjust = "left") +
      annotate(geom = "line", x = c(.43, 0.49), y = .11, colour = "#377EB8", size = 1.2) +
      annotate("text", x = 0.5, y = .11, size = 4.35, label = AUC_anno_2, color = "black", hjust = "left") +
      annotate(geom = "line", x = c(.43, 0.49), y = .05, colour = "#4DAF4A", size = 1.2) +
      annotate("text", x = 0.5, y = .05, size = 4.35, label = AUC_anno_3, color = "black", hjust = "left") +
      labs(x = "1-Specificity", y = "Sensitivity") +
      theme(
        axis.text.x = element_text(face = "plain", size = 12, color = "black"),
        axis.text.y = element_text(face = "plain", size = 12, color = "black"),
        axis.title.x = element_text(face = "plain", size = 14, color = "black"),
        axis.title.y = element_text(face = "plain", size = 14, color = "black")
      )

    if (do_ROC_CI) {
      ROC_95CI <- c()
      for (x in 1:3) {
        low <- round(ROC_data$AUC[[x]] - ROC_data[["inference"]][["vect_sd_1"]][x], digits = 3)
        up <- round(ROC_data$AUC[[x]] + ROC_data[["inference"]][["vect_sd_1"]][x], digits = 3)

        ROC_95CI[x] <- str_glue("(95%CI = {low} - {up})")
      }

      ROC_1 <- ggplot(data = rocFORplot) +
        geom_line(aes(x = year1x, y = year1y), size = 1.2, color = "#E83D3F") +
        geom_line(aes(x = year2x, y = year2y), size = 1.2, color = "#377EB8") +
        geom_line(aes(x = year3x, y = year3y), size = 1.2, color = "#4DAF4A") +
        geom_abline(slope = 1, intercept = 0, color = "grey", size = 1, linetype = 2) +
        egg::theme_article() +
        theme(plot.background = element_rect(fill = "white")) +
        annotate(geom = "line", x = c(.43, 0.49), y = .03 + 0.05 * 5, colour = "#E83D3F", size = 1.2) +
        annotate("text", x = 0.5, y = .03 + 0.05 * 5, size = 4.35, label = AUC_anno_1, color = "black", hjust = "left") +
        annotate("text", x = 0.5, y = .03 + 0.05 * 4, size = 4.35, label = ROC_95CI[1], color = "black", hjust = "left") +
        annotate(geom = "line", x = c(.43, 0.49), y = .03 + 0.05 * 3, colour = "#377EB8", size = 1.2) +
        annotate("text", x = 0.5, y = .03 + 0.05 * 3, size = 4.35, label = AUC_anno_2, color = "black", hjust = "left") +
        annotate("text", x = 0.5, y = .03 + 0.05 * 2, size = 4.35, label = ROC_95CI[2], color = "black", hjust = "left") +
        annotate(geom = "line", x = c(.43, 0.49), y = .03 + 0.05, colour = "#4DAF4A", size = 1.2) +
        annotate("text", x = 0.5, y = .03 + 0.05, size = 4.35, label = AUC_anno_3, color = "black", hjust = "left") +
        annotate("text", x = 0.5, y = .03, size = 4.35, label = ROC_95CI[3], color = "black", hjust = "left") +
        labs(x = "1-Specificity", y = "Sensitivity") +
        theme(
          axis.text.x = element_text(face = "plain", size = 12, color = "black"),
          axis.text.y = element_text(face = "plain", size = 12, color = "black"),
          axis.title.x = element_text(face = "plain", size = 14, color = "black"),
          axis.title.y = element_text(face = "plain", size = 14, color = "black")
        )
    }

    res <- list(NA, NA, NA, NA)
    res[[1]] <- KM
    res[[2]] <- ROC_1
    res[[3]] <- score_res

    p_auc_res <- tibble(
      "KM_p" = logrank_pvalue,
      "auc_1" = ROC_data$AUC[[1]], "auc_2" = ROC_data$AUC[[2]], "auc_3" = ROC_data$AUC[[3]],
      "HR" = round(HR, 2)
    )

    p <- ggarrange(res[[1]]$plot + theme(axis.title.x = element_blank()), res[[1]]$table, ncol = 1, align = "v", heights = c(0.7, 0.3))
    # print(cowplot::plot_grid(p, res[[2]], rel_widths = c(0.95, 1.05), scale = c(1, .97)))

    res[[4]] <- p_auc_res
  } else {
    res <- list(NA, NA, NA, NA)
  }

  return(res)
}

#' @TODO # 模型构建
#' @title  # 模型构建
#' @description
#'  - 包括对单因素cox回归分析，与lasso建模。
#'  - 调用子函数`univariate_Cox`,`lasso`,`KM_ROC_curve`.
#' @param exp 表达谱数据
#' @param clinical 临床信息文件
#' @param pval_univ_cox 单因素回归卡p值
#' @param seed 种子设置
#' @param saveplot 是否保存图片
#' @param output_dir 图片结果、lasso结果输出目录
#' @param surtime_unit 生存数据时间单位，12或者365
#' @param var_name 如果是NULL那么使用exp变量名字命名文件与文件夹
#' @export
#' @return  a list
#'  - 元素1：单因素回归结果，*data.frame*；
#'  - 元素2：lasso回归结果，*list*，四个元素：
#'    - 包含回归系数的data.frame
#'    - fit，可用于plot(fit)
#'    - cvfit，可用于plot(cvfit)
#'    - lasso回归系数图
#'  - 元素3：训练集结果*list*
#'    - KM曲线,gg对象
#'    - ROC曲线,gg对象
#'    - 训练集结果data.frame
#'
#' @author *WYK*
univCox2KM_ROC <- function(exp = NULL, clinical = NULL, pval_univ_cox = 0.1,
                           seed = 1110, saveplot = F, output_dir = "./",
                           surtime_unit = c(12, 365), ROC_time_break = c(1, 3, 5),
                           var_name = NULL, savetiff = F, best_cut = F,
                           by_group = T, do_ROC_CI = F) {
  if (is.null(var_name)) {
    var_name <- deparse(substitute(exp))
  }

  univcox_res <- univariate_Cox(exp = exp, clinical = clinical, pval_univ_cox = pval_univ_cox, by_group = by_group)
  print(sprintf("There're %s genes after univ_cox.", nrow(univcox_res)))

  if (nrow(univcox_res) > 1) {
    lasso_res <- lasso(
      clinical = clinical,
      exp = exp,
      unicox_res_gene = univcox_res$gene,
      seed = seed,
      saveplot = saveplot,
      output_dir = output_dir
    ) # lasso
    if (nrow(lasso_res[[1]]) > 1) {
      training <- KM_ROC_curve(
        model_coef = lasso_res[[1]],
        exp = exp,
        clinical = clinical,
        surtime_unit = surtime_unit,
        saveplot = saveplot,
        output_dir = output_dir,
        var_name = var_name,
        ROC_time_break = ROC_time_break,
        do_ROC_CI = do_ROC_CI,
        savetiff = savetiff,
        best_cut = best_cut
      )

      tmp <- list(univcox_res = univcox_res, lasso_res = lasso_res, training = training)
    } else {
      tmp <- list("univcox_res" = univcox_res, "lasso_res" = NA, "training" = NA)
    }
  } else {
    tmp <- list("univcox_res" = NA, "lasso_res" = NA, "training" = NA)
  }

  return(tmp)
}



#' @TODO 不同基因集差异分析
#' @title  # 不同基因集差异分析
#' @description  调用多个子函数
#' 先对表达谱进行差异分析，然后与兴趣基因集取交集，得到差异表达的兴趣基因，后续单因素cox回归 + lasso-cox回归建模.
#'
#' @param exp 表达谱数据,是一个变量名
#' @param clinical 临床信息表数据，要求 sample status time做列名。
#' @details 临床数据的样本名，需要和表达谱中样本名一致，或者绝大部分都相同
#' @param intrest_gene 聚类使用基因集
#' @param genelist 兴趣基因集，字符串向量
#' @param log2FC_value 差异基因阈值，绝对值
#' @param DEG_adjusted_pvalue 差异基因阈值 校正后p值
#' @param DEG_pvalue 差异基因阈值，p值 ，如果不为0，优先使用p值,如果是`NULL`，则使用校正后p值
#' @param pval_univ_cox 单因素回归P值
#' @param output_dir 结果输出路径，会在该路径下生成一个output文件夹，存放所有结果，需要以`"/"`结尾
#' @param seed 设置种子,默认1110
#' @param saveplot 是否保存图片，lasso结果图 与 训练集KM、ROC
#' @param surtime_unit 临床生存信息所用单位12或365
#' @param var_name 用于文件以及文件夹命名，需要是一个字符串。如果为NULL则使用随机字符做命名。
#' @param ROC_time_break ROC曲线时间断点默认1,3,5年。例如：`ROC_time_break = c(1,3,5)`
#' @param by_group 单因素cox回归分析是否使用二分组进行分析
#' @export
#'
#' @return   *list*
#'  - 元素1：*data.frame*，`deg_res`，阈值、genelist筛选后的差异基因结果
#'  - 元素2：*data.frame*，`deg_all_res`，差异分析全部结果
#'  - 元素3：*data.frame*，`univcox_res`，单因素cox回归分析结果
#'  - 元素4：*list*，`lasso_res`，子函数lasso分析结果
#'  - 元素5：*list*，`training`，训练集结果，包含km roc曲线以及训练集分组结果data.frame
#'
#' @author *WYK*
#'
#'
major_test1 <- function(exp = NULL, clinical = NULL,
                        intrest_gene = NULL,
                        group_infor = NULL, log2FC_value = 1,
                        DEG_adjusted_pvalue = 0.05, DEG_pvalue = NULL,
                        pval_univ_cox = 0.05, by_group = F,
                        surtime_unit = c(12, 365), ROC_time_break = c(1, 3, 5),
                        seed = 1110, saveplot = F,
                        output_dir = "./", var_name = NULL, do_ROC_CI = F) {
  if (length(var_name) == 0) {
    var_name <- paste0(sample(c(letters, LETTERS), 4), collapse = "")
  }

  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive = T) # 层级创建目录 类似mkdir -p
  } else {
    message(sprintf("'%s' is ready.", output_dir))
  }

  degs_list <- DEGs(exp = exp, group_infor = group_infor, output_dir = output_dir, var_name = var_name) # 计算差异基因

  deg_res <- model_data(
    degs = degs_list$degs,
    log2FC_value = log2FC_value,
    DEG_adjusted_pvalue = DEG_adjusted_pvalue,
    DEG_pvalue = DEG_pvalue,
    exp_1 = degs_list$exp_1,
    intrest_gene = intrest_gene
  ) # 按阈值筛选差异基因

  deg_gene_res_chara <- deg_res[[1]] %>%
    dplyr::select(gene) %>%
    pull() %>%
    unique()

  if (length(deg_gene_res_chara) > 1) {
    major_res <- univCox2KM_ROC(
      exp = exp[deg_gene_res_chara, ],
      clinical = clinical,
      pval_univ_cox = pval_univ_cox,
      seed = seed,
      saveplot = saveplot,
      output_dir = output_dir,
      surtime_unit = surtime_unit,
      var_name = var_name,
      by_group = by_group,
      ROC_time_break = ROC_time_break,
      do_ROC_CI = do_ROC_CI
    )

    if (length(major_res) > 1) {
      tmp3 <- list(
        deg_res = deg_res,
        deg_all_res = degs_list[[1]],
        univcox_res = major_res$univcox_res,
        lasso_res = major_res$lasso_res,
        training = major_res$training
      )
    } else {
      tmp3 <- list()
    }
  } else {
    tmp3 <- list()
  }

  return(tmp3)
}



#' @TODO 多因素cox回归计算HR
#' @title 多因素gene分析
#' @param exp 表达谱数据
#' @param model_coef_gene 基因
#' @param clinical 临床信息数据
#' @param output_dir 结果输出路径，需要以/结尾
#' @param saveplot 是否保存图片以及结果文件
#' @export
#' @return  a list
#' @author *WYK*
multi_Cox_gene <- function(exp = NULL, model_coef_gene = NULL, clinical = NULL, saveplot = F, output_dir = NULL) {
  exp <- exp[model_coef_gene, ]

  library(tidyverse)
  library(survival)
  # exp is only used to filter samples
  exp <- as.data.frame(t(exp))
  clinical <- clinical[na.omit(match(rownames(exp), clinical[, "sample"])), ]
  exp <- exp[na.omit(match(clinical[, "sample"], rownames(exp))), ]
  multicox_formulas <-
    as.formula(paste(
      "Surv(clinical$time, clinical$status)~",
      paste0(sep = "`", colnames(exp)[1:ncol(exp)], sep = "`", collapse = "+")
    ))
  # top 4 colunm：sample,time,status,riskgroup

  result_0 <- data.frame()

  result_0 <- coxph(
    formula = multicox_formulas,
    data = exp
  )

  result <- summary(result_0)

  result2 <- data.frame(
    var = rownames(result$coefficients),
    coef = result$coefficients[, 1],
    Hazard_Ratio = result$coefficients[, 2],
    p.value = result$coefficients[, 5],
    lower_.95 = result$conf.int[, "lower .95"],
    upper_.95 = result$conf.int[, "upper .95"]
  )

  result2 <- result2 %>%
    mutate(HR = paste0(
      round(Hazard_Ratio, 2),
      "(",
      round(lower_.95, 2),
      "-",
      round(upper_.95, 2),
      ")"
    ))
  str(result2)

  temp <- list()
  temp[[1]] <- result2
  temp[[2]] <- result_0


  p <- survminer::ggforest(
    model = temp[[2]],
    data = clinical,
    fontsize = 1,
    main = "Multivariable Analysis"
  )
  # print(p)

  temp[[3]] <- p

  if (saveplot) {
    if (!dir.exists(paste0(output_dir, "/multi_Cox_gene_forestplot/"))) {
      dir.create(paste0(output_dir, "/multi_Cox_gene_forestplot/"),
        mode = "7777", showWarnings = F, recursive = T
      )
    } else {
      (
        print("Dir is ready.")
      )
    }

    dir_now <- paste0(output_dir, "/multi_Cox_gene_forestplot/")
    cowplot::ggsave2(paste0(dir_now, "/multi_Cox_forestplot.pdf"),
      plot = p, width = 9, height = nrow(result2) * .7
    )
    # cowplot::ggsave2(paste0(dir_now,"/multi_Cox_forestplot.tiff"),
    #                  plot = p, width = 9, height = nrow(temp[[1]]) *.32, dpi = 300
    # )
    # cowplot::ggsave2(paste0(dir_now,"/multi_Cox_forestplot_dpi72.tiff"),
    #                  plot = p, width = 9, height = nrow(temp[[1]]) *.32, dpi = 72
    # )

    write.csv(x = result2, file = paste0(dir_now, "/multi_Cox_forestplot_gene.csv"), quote = F)
  } else {
    print("Save plot by Parameter 'saveplot = T'.")
  }

  temp[[3]] <- p

  names(temp) <- c("MUltiCox_res_df", "MUltiCox_res", "MultiCox_plot")

  return(temp)
}